function [trigger,triggerParams,triggerAdjusted]=ecogTriggersFindTriggersFromFile(analogueChannels2Use,ecogChannel,triggerValueRange,triggerDuration,interval2Analyze,checkResultVisually,pre,dur,scaleFac) %Find trigger events in a time series.
%[trigger,triggerParams,triggerAdjusted]=ecogTriggersFindTriggersFromFile(analogueChannels2Use,ecogChannel,triggerValueRange,triggerDuration,interval2Analyze,checkResultVisually,pre,dur) %Find trigger events in a time series.
%
%PURPOSE: Find trigger events in a time series loaded from file (e.g.
%created by ecogTDTData2Matlab) and return time series with trigger events
%marked.
%
%INPUT:
%analogueChannels2Use:  Files holding trigger timeseries (cell array of
%                       strings)
%ecogChannel:           File holding an ecog channel (string)
%triggerValueRange:     Range of values for the trigger event e.g. [3 4].
%                       One row per channel
%triggerDurationSec:    Max duration of a trigger in seconds e.g. 0.05.
%                       One value per channel
%                       This acts like a dead time. Only the first value
%                       in the range will be counted. If, however, the is
%                       actually longer than this interval an new trigger
%                       will be found at triggerduration+1
%intervalStartEndSec:   The intervall processed (in seconds)
%checkResultVisually:   Shows he timeseries around a trigger event.
%                       Sampling rate of trigger time series is adjusted to
%                       ecog time series
%pre:                   Time in seconds to show preceeding the trigger
%dur:                   Duration of interval to show around trigger.
%scaleFac:              A scalin factor to visualize the trigger timeseries
%
%OUTPUT:
%trigger:               A matrix with triggers marked. All other entries
%                       are zero. Rows correspond to timeseries. Original
%                       sampling rate is retained.
%triggerParams:         parameters of the original trigger timeseries
%triggerAdjusted:       A vector with triggers marked. Sampling rate is
%                       adjusted to ECoG-waves.
% EXAMPLE:
%[trigger,triggerParams,triggerAdjusted]=ecogTriggersFindTriggersFromFile({'ANIN3','ANIN4'},'Wav11',[3 4;3 4],[0.4 0.4],[4 188],1,0.002,0.004,1);

% 20110810 JR wrote it

%Load the analogue channel
for k=1:length(analogueChannels2Use)
    if k==1
        %load an ecog timeseries to fetch parameters
        load(ecogChannel)
        nSampInWav=length(data);
        waveParams=params;
        triggerAdjusted=zeros(length(analogueChannels2Use),length(data)); %original sampling rate trigger
    end
    % load the analogue timeseries
    load(analogueChannels2Use{k});
    triggerParams=params;
    %initialize some values and variables
    if k==1
        trigger=zeros(length(analogueChannels2Use),length(data)); %original sampling rate trigger
        intervalStartEndSample=nearly(interval2Analyze,(1:length(data))/triggerParams.sampFreq);
        minDurOfStimulusInSamples=round(triggerDuration(k)*triggerParams.sampFreq);
    end
    % set everything outside of the intervall of interest to zero and find
    % trigger events data is a vectr loaded from file. 
    data(1:intervalStartEndSample(1))=0;
    data(intervalStartEndSample(2):end)=0;
    threshData=double(data>triggerValueRange(k,1) & data<triggerValueRange(k,2));
    count=1;
    while count<intervalStartEndSample(2)
        if threshData(count)==1
            trigger(k,count)=1; %set a trigger
            count=count+minDurOfStimulusInSamples; % adavance to end of dead interval
        else
            count=count+1;
        end
    end
    
    %Check the triggers on full sampling rate. A remainder from debugging.
    %Not sure it it works
    if 0
        if checkResultVisually
            %check if all onsets are correct
            idx=find(trigger(k,:));
            figure;
            r=ceil(-pre*triggerParams.sampFreq:(-pre+dur)*triggerParams.sampFreq);
            timebase=r*1000/triggerParams.sampFreq;
            for m=1:length(idx)
                plot(timebase,[threshData(r+idx(m));data(r+idx(m));trigger(k,r+idx(m))*2]');
                input(num2str(m));
            end
        end
    end
        
    %Adjust trigger channel sampling rate to ecog sampling rate
    decFac=triggerParams.sampFreq/waveParams.sampFreq;
    oldIdx=find(trigger(k,:)>0);
    newIdx=round(oldIdx/decFac);
    for m=1:length(newIdx)
        triggerAdjusted(k,newIdx(m))=trigger(k,oldIdx(m));
    end
    
    if checkResultVisually
        
        % downsample the original trigger timeseries
        [B,A] = butter(3,1/(2*decFac),'low'); %factor to nyquits of downsampled sequence
        data=filtfilt(B,A,double(data));
        aninDown=data(1:decFac:end);
        
        idx=find(triggerAdjusted(k,:));
        figure;
        r=ceil(-pre*waveParams.sampFreq:(-pre+dur)*waveParams.sampFreq);
        timebase=r*1000/waveParams.sampFreq;
        for m=1:length(idx)
            plot(timebase,[aninDown(r+idx(m));triggerAdjusted(k,r+idx(m))*scaleFac]');
            title(['Channel:' analogueChannels2Use{k} ' trial # : ' num2str(m)])
            input(num2str(m));
        end
        if length(analogueChannels2Use)>k
            flag = 0;
            while flag==0
                r=input('Continue processing next set? [y]/n ','s');
                if strcmpi('y',r)
                    flag = 1;
                    % Next channel
                elseif strcmpi('n',r)
                    flag=0;
                    error(['User aborted. Current channel: ' analogueChannels2Use{k}]);
                else
                    disp('Inavlid answer. Try again.')
                    flag=0;
                end
            end
        end     
    end
    
end